Return to Main Page

\(\color{darkblue}{\textbf{sf}}\)

---

Used for processing spatial data and basic visualization

library(tidyverse)
library(sf)
library(raster)

---

\(\color{dodgerblue}{\textbf{Importing}}\)

\(\color{skyblue}{\textrm{- Reading in Data}}\)

Vector data (in any spatial form) can be read in with read_sf…

url <- "https://github.com/social-lorax/new_guides/raw/main/Spatial/Borough%20Boundaries.geojson"
boros <- read_sf(url) 

…or transformed using st_as_sf()

spatial <- read_csv("pathway/to/spatial.csv") %>% 
  st_as_sf()

Raster data is read in using raster() if it is single-band…

url <- "https://github.com/social-lorax/new_guides/blob/main/Spatial/canopy.tif?raw=true"
canopy <- raster(url)

…and brick() if it is multi-band

url <- "https://github.com/social-lorax/new_guides/blob/main/Spatial/manhattan.tif?raw=true"
manhattan <- brick(url)

---

\(\color{dodgerblue}{\textbf{Removing Geometry}}\)

Having a geometry can be annoying for other functions, so st_set_geometry(NULL) removes it

non_spatial <- spatial %>% 
  st_set_geometry(NULL)

---

\(\color{dodgerblue}{\textbf{Geocoding}}\)

boro_halls <- tibble(hall = c("New York City Hall", "Bronx County Courthouse", "Brooklyn Borough Hall",
                              "Queens Borough Hall", "Staten Island Borough Hall"),
                     street = c("New York City Hall, New York, NY 10007", 
                                "851 Grand Concourse, Bronx, NY 10451", 
                                "209 Joralemon St, Brooklyn, NY 11201",
                                "120-55 Queens Blvd, Queens, NY 11424", 
                                "10 Richmond Terrace, Staten Island, NY 10301"))

boro_halls

hall	street
New York City Hall	New York City Hall, New York, NY 10007
Bronx County Courthouse	851 Grand Concourse, Bronx, NY 10451
Brooklyn Borough Hall	209 Joralemon St, Brooklyn, NY 11201
Queens Borough Hall	120-55 Queens Blvd, Queens, NY 11424
Staten Island Borough Hall	10 Richmond Terrace, Staten Island, NY 10301

\(\color{skyblue}{\textrm{- Cleaning Messy Addresses}}\)

library(postmastr) #remotes::install_github("slu-openGIS/postmastr")

#Set-up
boro_halls <- boro_halls %>% 
  pm_identify(var = "street")

boro_halls_unique <- boro_halls %>% 
  pm_prep(var = "street", type = "street")

ZIP Code

#Check
boro_halls_unique %>% pm_postal_none()

pm.uid	pm.address

#Parse
boro_halls_unique <- boro_halls_unique %>% 
  pm_postal_parse()

State

#Check
nyDict <- pm_dictionary(locale = "us", type = "state", filter = "NY", case = c("title", "upper"))

boro_halls_unique %>% pm_state_none(dictionary = nyDict)

pm.uid	pm.address	pm.zip

#Parse
boro_halls_unique <- boro_halls_unique %>% 
  pm_state_parse()

City

#Check
defaultDict <- pm_dictionary(locale = "us", type = "city", filter = "NY")

specificDict <- pm_append(type = "city",
                          input = c("New York", "New York City", "NYC"     , "Bronx", "The Bronx",
                                    "Brooklyn", "Queens", "Jamaica", "Staten Island"),
                          output = c(NA       , "New York"     , "New York", NA     , "Bronx"    , 
                                     NA       , NA      , "Queens" , NA))

boro_halls_unique %>% pm_city_none(dictionary = specificDict)

pm.uid	pm.address	pm.state	pm.zip

#Parse
boro_halls_unique <- boro_halls_unique %>% 
  pm_city_parse(dictionary = specificDict)

House Number

#Check
boro_halls_unique %>% pm_house_none()

pm.uid	pm.address	pm.city	pm.state	pm.zip
1	New York City Hall	New York	NY	10007

#Parse after fixing or ignoring
boro_halls_unique <- boro_halls_unique %>% 
  pm_house_parse()

Street Direction (e.g., East Broadway)

#Check 
boro_halls_unique %>% pm_streetDir_none()

pm.uid	pm.address	pm.house	pm.city	pm.state	pm.zip
1	New York City Hall	NA	New York	NY	10007
2	Grand Concourse	851	Bronx	NY	10451
3	Joralemon St	209	Brooklyn	NY	11201
4	Queens Blvd	120-55	Queens	NY	11424
5	Richmond Terrace	10	Staten Island	NY	10301

Moving on since none have a prefix

Street Suffix (e.g., St., Ave.)

#Check
boro_halls_unique %>% pm_streetSuf_none()

pm.uid	pm.address	pm.house	pm.city	pm.state	pm.zip
1	New York City Hall	NA	New York	NY	10007
2	Grand Concourse	851	Bronx	NY	10451

#Parse after fixing or ignoring
boro_halls_unique <- boro_halls_unique %>% 
  pm_streetSuf_parse()

Street

#Street is assumed 
boro_halls_unique <- boro_halls_unique %>% 
  pm_street_parse()

Reassemble

boro_halls_parsed <- pm_replace(boro_halls_unique, source = boro_halls) %>% 
  pm_rebuild(output = "full", keep_parsed = "no")

boro_halls_parsed

hall	street	pm.address
New York City Hall	New York City Hall, New York, NY 10007	New York City Hall New York NY 10007
Bronx County Courthouse	851 Grand Concourse, Bronx, NY 10451	851 Grand Concourse Bronx NY 10451
Brooklyn Borough Hall	209 Joralemon St, Brooklyn, NY 11201	209 Joralemon St Brooklyn NY 11201
Queens Borough Hall	120-55 Queens Blvd, Queens, NY 11424	120-55 Queens Blvd Queens NY 11424
Staten Island Borough Hall	10 Richmond Terrace, Staten Island, NY 10301	10 Richmond Ter Staten Island NY 10301

\(\color{skyblue}{\textrm{- Geolocating}}\)

library(tidygeocoder)

processed <- boro_halls_parsed %>% 
  geocode(address = pm.address, method = "osm", verbose = TRUE, unique_only = TRUE)

boro_halls <- boro_halls_parsed %>% 
  full_join(processed, by = c("pm.address" = "address" )) %>% 
  st_as_sf(coords = c("long", "lat"), crs = 4326)

boro_halls

hall	street	pm.address	geometry
New York City Hall	New York City Hall, New York, NY 10007	New York City Hall New York NY 10007	POINT (-74.00594 40.71274)
Bronx County Courthouse	851 Grand Concourse, Bronx, NY 10451	851 Grand Concourse Bronx NY 10451	POINT (-73.92332 40.82611)
Brooklyn Borough Hall	209 Joralemon St, Brooklyn, NY 11201	209 Joralemon St Brooklyn NY 11201	POINT (-73.9903 40.69276)
Queens Borough Hall	120-55 Queens Blvd, Queens, NY 11424	120-55 Queens Blvd Queens NY 11424	POINT (-73.8283 40.71378)
Staten Island Borough Hall	10 Richmond Terrace, Staten Island, NY 10301	10 Richmond Ter Staten Island NY 10301	POINT (-74.07609 40.64242)

\(\color{skyblue}{\textrm{- NYC Geocoding}}\)

#install.packages("remotes")
remotes::install_github("austensen/geoclient")
geoclient::geoclient_api_key("key")

geoclient::geo_search_data(boro_halls_parsed, street) %>% 
  transmute(address = input_location, 
            no_results, 
            bbl,
            bin = buildingIdentificationNumber,
            tract = censusTract2010,
            block = censusBlock2010,
            cd_name = communityDistrict,
            cd_num = communityDistrictNumber,
            puma = pumaCode,
            nta_name = ntaName,
            nta_num = nta,
            citycouncil = cityCouncilDistrict,
            assembly = assemblyDistrict,
            congressional = congressionalDistrict,
            latitude,
            longitude)

address	no_results	bbl	bin	tract	block	cd_name	cd_num	puma	nta_name	nta_num	citycouncil	assembly	congressional	latitude	longitude
New York City Hall, New York, NY 10007	FALSE	1001220001	1079147	31	1021	101	01	03810	SoHo-TriBeCa-Civic Center-Little Italy	MN24	01	66	10	40.71281	-74.00610
851 Grand Concourse, Bronx, NY 10451	FALSE	2024680001	2002869	63	3000	204	04	03708	West Concourse	BX63	08	84	15	40.82627	-73.92297
209 Joralemon St, Brooklyn, NY 11201	FALSE	3001390001	3000256	11	1003	302	02	04004	DUMBO-Vinegar Hill-Downtown Brooklyn-Boerum Hill	BK38	33	52	07	40.69248	-73.99049
120-55 Queens Blvd, Queens, NY 11424	FALSE	4022740002	4052812	216	1008	409	09	04111	Kew Gardens	QN60	29	27	06	40.71332	-73.82858
10 Richmond Terrace, Staten Island, NY 10301	FALSE	5000070001	5000063	3	1004	501	01	03903	West New Brighton-New Brighton-St. George	SI22	49	61	11	40.64244	-74.07533

---

\(\color{dodgerblue}{\textbf{Projection}}\)

\(\color{skyblue}{\textrm{- Checking}}\)

For vector data, use st_crs()

boro_halls %>% st_crs()

Coordinate Reference System:
  User input: EPSG:4326 
  wkt:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    USAGE[
        SCOPE["Horizontal component of 3D system."],
        AREA["World."],
        BBOX[-90,-180,90,180]],
    ID["EPSG",4326]]

For raster data, just use crs()

canopy %>% crs()

CRS arguments:
 +proj=aea +lat_0=23 +lon_0=-96 +lat_1=29.5 +lat_2=45.5 +x_0=0 +y_0=0
+ellps=GRS80 +units=m +no_defs 

\(\color{skyblue}{\textrm{- Reprojecting}}\)

Name	Code
World Geodetic System 1984 (WGS 84)	4326
NAD83 / New York Long Island	2263

For vector data, use st_transform()

boro_halls_2263 <- boro_halls %>% 
  st_transform(2263)

boros_2263 <- boros %>% 
  st_transform(2263)

For raster data, use projectRaster()

canopy_2263 <- canopy %>% 
  projectRaster(crs = 2263)

---

\(\color{dodgerblue}{\textbf{Simplifying}}\)

For vector data, use st_simplify() (you need to transform out of lat/long first but can transform back after)

object.size(boros)

1372240 bytes

boros_simple <- boros %>% 
  st_transform(2263) %>% 
  st_simplify(preserveTopology = TRUE,
              dTolerance = 100) %>% 
  st_transform(4326)

object.size(boros_simple)

116320 bytes

str_c(round(100 * (116320 - 1372240) / 1372240, 2), "%")

[1] "-91.52%"

plot(boros$geometry)

plot(boros_simple$geometry)

For raster data, use aggregate()

manhattan_simple <- manhattan %>% 
  aggregate(fact = 5, fun = mean)

res(manhattan)

[1] 29.98979 30.00062

res(manhattan_simple)

[1] 149.9489 150.0031

ncell(manhattan)

[1] 619173

ncell(manhattan_simple)

[1] 24955

plotRGB(manhattan)

plotRGB(manhattan_simple)

---

\(\color{dodgerblue}{\textbf{Plotting}}\)

\(\color{skyblue}{\textrm{- Basic}}\)

ggplot() + 
  geom_sf(data = boros, aes(fill = population)) + 
  geom_sf(data = boro_halls, color = "red", size = 2) + 
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        rect = element_blank()) + 
  scale_fill_viridis_c()

\(\color{skyblue}{\textrm{- Basemap}}\)

library(ggspatial)

Options:

Bounding box to zoom:

boro_halls %>% 
  ggplot() + 
  annotation_map_tile(type = "cartolight", zoomin = 0) +
  geom_sf(color = "red", size = 1.5) + 
  coord_sf(xlim = c(-74.2591, -73.7002), 
           ylim = c(40.4774, 40.9162), 
           expand = FALSE) +
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        rect = element_blank())

---

\(\color{dodgerblue}{\textbf{Analysis}}\)

\(\color{skyblue}{\textrm{- Centroids}}\)

st_centroid() creates a point at the center of each polygon

boros_centers <- boros_2263 %>% 
  st_centroid()

ggplot() + 
  geom_sf(data = boros_2263) +
  geom_sf(data = boros_centers, color = "blue", size = 2) + 
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        rect = element_blank()) 

\(\color{skyblue}{\textrm{- Buffering}}\)

First check the projection to get the units you are working with (here, US survey foot, which is about 0.999998 feet)

boro_halls_2263 %>% st_crs()

Coordinate Reference System:
  User input: EPSG:2263 
  wkt:
PROJCRS["NAD83 / New York Long Island (ftUS)",
    BASEGEOGCRS["NAD83",
        DATUM["North American Datum 1983",
            ELLIPSOID["GRS 1980",6378137,298.257222101,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4269]],
    CONVERSION["SPCS83 New York Long Island zone (US Survey feet)",
        METHOD["Lambert Conic Conformal (2SP)",
            ID["EPSG",9802]],
        PARAMETER["Latitude of false origin",40.1666666666667,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8821]],
        PARAMETER["Longitude of false origin",-74,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8822]],
        PARAMETER["Latitude of 1st standard parallel",41.0333333333333,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8823]],
        PARAMETER["Latitude of 2nd standard parallel",40.6666666666667,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8824]],
        PARAMETER["Easting at false origin",984250,
            LENGTHUNIT["US survey foot",0.304800609601219],
            ID["EPSG",8826]],
        PARAMETER["Northing at false origin",0,
            LENGTHUNIT["US survey foot",0.304800609601219],
            ID["EPSG",8827]]],
    CS[Cartesian,2],
        AXIS["easting (X)",east,
            ORDER[1],
            LENGTHUNIT["US survey foot",0.304800609601219]],
        AXIS["northing (Y)",north,
            ORDER[2],
            LENGTHUNIT["US survey foot",0.304800609601219]],
    USAGE[
        SCOPE["Engineering survey, topographic mapping."],
        AREA["United States (USA) - New York - counties of Bronx; Kings; Nassau; New York; Queens; Richmond; Suffolk."],
        BBOX[40.47,-74.26,41.3,-71.8]],
    ID["EPSG",2263]]

st_buffer() will buffer around a feature (though you cannot add units, so you have to know!)

buffers <- boro_halls_2263 %>% 
  st_buffer(dist = 5280)

ggplot() + 
  geom_sf(data = boros_2263) + 
  geom_sf(data = buffers, fill = "red", alpha = 0.5) + 
  geom_sf(data = boro_halls, color = "black", size = 2) + 
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        rect = element_blank()) 

\(\color{skyblue}{\textrm{- Distances}}\)

st_distance() returns a matrix of the distance between each input to each output

boro_halls_2263 %>% 
  st_distance(boros_centers)

Units: [US_survey_foot]
         [,1]     [,2]      [,3]     [,4]      [,5]
[1,] 63944.00 25834.72  63119.49 29549.81 51982.247
[2,] 18452.27 21556.31 109797.41 66425.85 52004.377
[3,] 67591.47 31438.44  60926.31 21084.47 47928.670
[4,] 51683.39 44898.17 102397.11 41619.33  3501.901
[5,] 96115.36 57654.35  31051.68 35598.83 75287.273

More processing is needed to find, for example, the closest point

hall_center <- boro_halls_2263 %>% 
  st_distance(boros_centers) %>%
  as_tibble()

names(hall_center) <- boros_centers$boro_name

hall_center %>% 
  mutate(hall = boro_halls_2263$hall) %>% 
  gather(-hall, key = boro_center, value = distance) %>% 
  mutate(distance = as.double(distance)) %>% 
  arrange(hall, distance) %>% 
  distinct(hall, .keep_all = TRUE) %>% 
  full_join(boro_halls_2263, by = "hall")

hall	boro_center	distance	street	pm.address	geometry
Bronx County Courthouse	Bronx	18452.274	851 Grand Concourse, Bronx, NY 10451	851 Grand Concourse Bronx NY 10451	POINT (1005472 240261.9)
Brooklyn Borough Hall	Brooklyn	21084.472	209 Joralemon St, Brooklyn, NY 11201	209 Joralemon St Brooklyn NY 11201	POINT (986940.8 191666.1)
New York City Hall	Manhattan	25834.723	New York City Hall, New York, NY 10007	New York City Hall New York NY 10007	POINT (982603.5 198946.8)
Queens Borough Hall	Queens	3501.901	120-55 Queens Blvd, Queens, NY 11424	120-55 Queens Blvd Queens NY 11424	POINT (1031849 199373.6)
Staten Island Borough Hall	Staten Island	31051.678	10 Richmond Terrace, Staten Island, NY 10301	10 Richmond Ter Staten Island NY 10301	POINT (963132.8 173336.4)

\(\color{skyblue}{\textrm{- Area}}\)

st_area() returns the area of each polygon (make sure to check the CRS for the units!)

boros_2263 %>% 
  st_area() %>% 
  as_tibble() %>% 
  transmute(boro = boros_2263$boro_name,
            area = as.double(value)) 

boro	area
Bronx	1187195250
Manhattan	636575581
Staten Island	1623637770
Brooklyn	1934173752
Queens	3040206185

\(\color{skyblue}{\textrm{- Length}}\)

st_length() returns the length of a line

library(tigris)

roads <- tigris::roads(state = "NY", county = c("Bronx", "Kings", "New York", "Queens", "Richmond"), progress_bar = FALSE) %>% 
  st_transform(2263)

bedford <- roads %>% 
  filter(str_detect(FULLNAME, "Bedford Ave"),
         !(LINEARID %in% c("110425921352", "110444047995"))) %>% 
  st_union()

ggplot() + 
  geom_sf(data = boros_2263) + 
  geom_sf(data = bedford, color = "yellow", size = 1) + 
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        rect = element_blank())

bedford %>% 
  st_length()

53673.32 [US_survey_foot]

\(\color{skyblue}{\textrm{- Dissolving}}\)

st_union() takes multiple entities and creates a single entity

mn_buffer <- boro_halls_2263 %>% 
  filter(hall == "New York City Hall") %>% 
  st_buffer(dist = 5000) %>% 
  st_transform(4326)

bk_buffer <- boro_halls_2263 %>% 
  filter(hall == "Brooklyn Borough Hall") %>% 
  st_buffer(dist = 5000) %>% 
  st_transform(4326)

ggplot() + 
  annotation_map_tile(type = "cartolight", zoomin = 0) +
  geom_sf(data = mn_buffer, fill = "red", color = "black", alpha = 0.5, ) + 
  geom_sf(data = bk_buffer, fill = "blue", color = "black", alpha = 0.5) + 
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        rect = element_blank()) 

joined_buffer <- mn_buffer %>% 
  st_union(bk_buffer)

ggplot() + 
  annotation_map_tile(type = "cartolight", zoomin = 0) +
  geom_sf(data = joined_buffer, fill = "purple", color = "black", alpha = 0.5, ) + 
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        rect = element_blank()) 

\(\color{skyblue}{\textrm{- Spatial Joins}}\)

st_join() returns everything, there is also st_intersects() and st_contains for more control

rand_points <- tibble(point = c("A", "B", "C", "D"),
                      lat = c(40.8845, 40.7202, 40.5905, 40.6366),
                      long = c(-73.7477, -73.9846, -74.0572, -73.7812)) %>% 
  st_as_sf(coords  = c("long", "lat"), crs = 4326) %>% 
  st_transform(crs = 2263)

rand_points %>% 
  st_join(boros_2263)

point	geometry	boro_code	boro_name	shape_area	shape_leng
A	POINT (1054014 261625.6)	NA	NA	NA	NA
B	POINT (988518.8 201665.2)	1	Manhattan	636575908.328	359804.890798
C	POINT (968363.5 154416.7)	NA	NA	NA	NA
D	POINT (1044977 171282.8)	4	Queens	3040205398.17	900270.512975

ggplot() + 
  geom_sf(data = boros_2263) + 
  geom_sf(data = rand_points, color = "black", size = 2) + 
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        rect = element_blank()) 

bk <- boros_2263 %>% filter(boro_name == "Brooklyn")

bk_intersect <- bk %>% st_intersects(roads)
bk_intersect <- roads[bk_intersect[[1]],]

bk_contained <- bk %>% st_contains(roads)
bk_contained <- roads[bk_contained[[1]],]

ggplot() + 
  geom_sf(data = bk) + 
  geom_sf(data = bk_intersect, color = "blue", size = 1) + 
  geom_sf(data = bk_contained, color = "yellow", size = 1, alpha = 0.5) + 
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        rect = element_blank()) 

\(\color{skyblue}{\textrm{- Clipping}}\)

st_intersection() clips one entity based on how much of it intersects with another entity

boro_clips <- buffers %>% 
  filter(hall == "Bronx County Courthouse") %>% 
  st_intersection(boros_2263)

ggplot() + 
  annotation_map_tile(type = "cartolight", zoomin = 0) +
  geom_sf(data = boro_clips, fill = "purple", color = "black", alpha = 0.5, ) + 
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        rect = element_blank()) 

\(\color{skyblue}{\textrm{- Bounding Box}}\)

st_make_grid() creates the smallest box possible containing all elements

boro_box <- boro_halls %>% 
  st_make_grid(n = 1)

ggplot() + 
  geom_sf(data = boros) + 
  geom_sf(data = boro_box, fill = "orange", alpha = 0.5) + 
  geom_sf(data = boro_halls, color = "black", size = 2) + 
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        rect = element_blank()) 

\(\color{skyblue}{\textrm{- Convex Hull}}\)

st_convex_hull() creates the smallest shape possible containing all elements

hull <- boro_halls %>% 
  st_union() %>% 
  st_convex_hull()

ggplot() + 
  geom_sf(data = boros) + 
  geom_sf(data = hull, fill = "orange", alpha = 0.5) + 
  geom_sf(data = boro_halls, color = "black", size = 2) + 
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        rect = element_blank()) 

\(\color{skyblue}{\textrm{- Masking Rasters}}\)

mask() and crop() together clip a raster

canopy_2263 %>% 
  raster::plot()

canopy_2263 %>% 
  mask(boros_2263) %>% 
  crop(boros_2263) %>% 
  raster::plot()

\(\color{skyblue}{\textrm{- Overlaying Rasters}}\)

overlay() with a function can manipulate a raster

f <- function(rast) {
  rast > 60 & rast < 100
}

canopy_2263 %>% 
  overlay(fun = f) %>%
  raster::plot()

\(\color{skyblue}{\textrm{- Extracting Rasters}}\)

extract() pulls aggregated raster information

mean_canopy <- canopy_2263 %>% 
  extract(boros_2263, fun = mean)

median_canopy <- canopy_2263 %>% 
  extract(boros_2263, fun = median)

max_canopy <- canopy_2263 %>% 
  extract(boros_2263, fun = max)

min_canopy <- canopy_2263 %>% 
  extract(boros_2263, fun = min)

boros %>% 
  mutate(mean_canopy = mean_canopy,
         median_canopy = median_canopy,
         min_canopy = min_canopy,
         max_canopy = max_canopy)

boro_code	boro_name	shape_area	shape_leng	geometry	population	mean_canopy	median_canopy	min_canopy	max_canopy
2	Bronx	1187193513.84	463868.937681	MULTIPOLYGON (((-73.89681 4…	2405464	13.062493	4.7337441	0	84.98292
1	Manhattan	636575908.328	359804.890798	MULTIPOLYGON (((-74.01093 4…	2736074	7.252479	2.1085426	0	76.45837
5	Staten Island	1623635734.34	325929.794027	MULTIPOLYGON (((-74.05051 4…	1694251	21.885133	15.4100775	0	87.58646
3	Brooklyn	1934174208.45	728195.197649	MULTIPOLYGON (((-73.86327 4…	1472654	3.652571	0.9797716	0	65.93524
4	Queens	3040205398.17	900270.512975	MULTIPOLYGON (((-73.82645 4…	495747	6.911773	2.9477946	0	86.00687

\(\color{darkblue}{\textbf{Leaflet}}\)

---

Used for creating interactive maps

library(leaflet)
library(leafpop)
library(leaflet.mapboxgl) # remotes::install_github("rstudio/leaflet.mapboxgl")

---

\(\color{dodgerblue}{\textbf{Basic Map}}\)

leaflet() %>% 
  addTiles() %>%
  addMarkers(lng = -74.000060, lat = 40.730910)

---

\(\color{dodgerblue}{\textbf{Zoom}}\)

Far (zoom = 10)

leaflet() %>% 
  addTiles() %>%
  setView(lng = -74.000060, lat = 40.730910, zoom = 10) %>% 
  addMarkers(lng = -74.000060, lat = 40.730910)

Close (zoom = 100)

leaflet() %>% 
  addTiles() %>%
  setView(lng = -74.000060, lat = 40.730910, zoom = 100) %>% 
  addMarkers(lng = -74.000060, lat = 40.730910)

---

\(\color{dodgerblue}{\textbf{Markers}}\)

boro_hall_points <- boro_halls %>% 
  st_transform(4326) %>% 
  st_coordinates() 

boro_halls <- boro_halls %>%  
  st_set_geometry(NULL) %>% 
  bind_cols(tibble(long = boro_hall_points[,1],
                   lat = boro_hall_points[,2]))

\(\color{skyblue}{\textrm{- Pin}}\)

boro_halls %>% 
  leaflet() %>% 
  addTiles() %>%
  addMarkers(lng = ~long, lat = ~lat)

\(\color{skyblue}{\textrm{- Colored Circle}}\)

boro_halls %>% 
  leaflet() %>% 
  addTiles() %>%
  addCircleMarkers(radius = 6,
                   color = "red",
                   fillOpacity = 0.8,
                   stroke = FALSE)

---

\(\color{dodgerblue}{\textbf{Labels}}\)

ny <- paste(sep = "<br/>",
            "<b><a href='https://www1.nyc.gov/'>City Hall</a></b>",
                 "City Hall Park",
                 "New York, NY 10007")

bx <- paste(sep = "<br/>",
            "<b><a href='https://bronxboropres.nyc.gov/'>Bronx Borough Hall</a></b>",
                 "851 Grand Concourse",
                 "Bronx, NY 10451")

bk <- paste(sep = "<br/>",
            "<b><a href='https://www.brooklyn-usa.org/'>Brooklyn Borough Hall</a></b>",
                 "209 Joralemon St",
                 "Brooklyn, NY 11201")

qn <- paste(sep = "<br/>",
            "<b><a href='https://queensbp.org/'>Queens Borough Hall</a></b>",
                 "120-55 Queens Blvd",
                 "Queens, NY 11424")

si <- paste(sep = "<br/>",
            "<b><a href='https://www.statenislandusa.com/'>Staten Island Borough Hall</a></b>",
                 "10 Richmond Terrace",
                 "Staten Island, NY 10301")

boro_halls <- boro_halls %>% 
  mutate(label = c(ny, bx, bk, qn, si))

\(\color{skyblue}{\textrm{- Popups}}\)

boro_halls %>% 
  leaflet() %>% 
  addTiles() %>%
  addMarkers(lng = ~long, lat = ~lat,
             popup = ~label)

\(\color{skyblue}{\textrm{- Hover}}\)

labels <- sprintf(str_glue("{boro_halls$hall}<br>{boro_halls$street}")) %>% 
  lapply(htmltools::HTML)

boro_halls %>% 
  leaflet() %>% 
  addTiles() %>%
  addMarkers(lng = ~long, lat = ~lat,
             label = labels)

\(\color{skyblue}{\textrm{- Formatting}}\)

leaflet() %>% 
  addTiles() %>% 
  setView(-118.456554, 34.09, 13) %>%
  addMarkers(lng = -118.456554, lat = 34.105,
             label = "Default Label",
             labelOptions = labelOptions(noHide = T)) %>%
  addMarkers(lng = -118.456554, lat = 34.095,
             label = "Label w/o surrounding box",
             labelOptions = labelOptions(noHide = T, 
                                         textOnly = TRUE)) %>%
  addMarkers(lng = -118.456554, lat = 34.085,
             label = "label w/ textsize 15px",
             labelOptions = labelOptions(noHide = T, 
                                         textsize = "15px")) %>%
  addMarkers(lng = -118.456554, lat = 34.075,
             label = "Label w/ custom CSS style",
             labelOptions = labelOptions(noHide = T, 
                                         direction = "bottom",
                                         style = list("color" = "red",
                                                      "font-family" = "serif",
                                                      "font-style" = "italic",
                                                      "box-shadow" = "3px 3px rgba(0,0,0,0.25)",
                                                      "font-size" = "12px",
                                                      "border-color" = "rgba(0,0,0,0.5)")))

---

\(\color{dodgerblue}{\textbf{Polygons}}\)

#First getting data from the census to plot (see section on census data)
library(tidycensus)
population <- get_acs(geography = "county",
                      state = "NY",
                      variables = c(population = "B01003_001"),
                      year = 2019, 
                      survey = "acs5",
                      geometry = TRUE,
                      progress_bar = FALSE) %>% 
  mutate(pretty_estimate = scales::comma(estimate),
         legend_value = estimate / 10000) %>% 
  st_as_sf()

pal <- colorBin("viridis", domain = population$legend_value)

labels <- sprintf("<strong>%s</strong><br/>Population: %s",
                  population$NAME, population$pretty_estimate) %>% 
  lapply(htmltools::HTML)

population %>% 
  leaflet() %>% 
  addTiles() %>% 
  addPolygons(fillColor = ~pal(legend_value),
              fillOpacity = 0.7,
              color = "white",
              weight = 1,
              opacity = 1,
              smoothFactor = 0.5,
              highlight = highlightOptions(weight = 5,
                                           color = "white",
                                           fillOpacity = 0.7,
                                           bringToFront = TRUE),
              label = labels,
              labelOptions = labelOptions(style = list("font-weight" = "normal", 
                                                       padding = "3px 8px"),
                                          textsize = "15px",
                                          direction = "auto")) %>% 
  addLegend(pal = pal, 
            values = ~legend_value, 
            opacity = 0.7, 
            title = "Population (10K)",
            position = "bottomright")

Return to Main Page